library(dplyr) 
library(org.Hs.eg.db) 
library(clusterProfiler) 
library(DOSE)
library(ggplot2) 
library(RColorBrewer)

diff <-rownames(C3) 


diff_entrez <- bitr(diff,
                    fromType = "SYMBOL",#现有的ID类型
                    toType = "ENTREZID",#需转换的ID类型
                    OrgDb = "org.Hs.eg.db")
head(diff_entrez) #提示少量无法映射，不同数据库间ID转换存在少量缺失是正常现象


#KEGG富集分析(超几何分布检验)：
go<-enrichGO(gene = diff_entrez$ENTREZID, OrgDb = "org.Hs.eg.db", ont="all",pvalueCutoff = 0.1)

KEGG_diff <- enrichKEGG(gene = diff_entrez$ENTREZID,
                        organism = "hsa", #物种Homo sapiens (智人)
                        pvalueCutoff = 0.05,
                        qvalueCutoff = 0.05,
                        pAdjustMethod = "BH",
                        minGSSize = 10,
                        maxGSSize = 500)

#将ENTREZ重转为symbol：
KEGG_diff <- setReadable(KEGG_diff,
                         OrgDb = org.Hs.eg.db,
                         keyType = "ENTREZID")
View(KEGG_diff@result)


#计算Rich Factor（富集因子）：
KEGG_diff2 <- mutate(KEGG_diff,
                     RichFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))

#计算Fold Enrichment（富集倍数）：
KEGG_diff2 <- mutate(KEGG_diff2, FoldEnrichment = parse_ratio(GeneRatio) / parse_ratio(BgRatio))
KEGG_diff2@result$RichFactor[1:6]
KEGG_diff2@result$FoldEnrichment[1:6]


library(enrichplot)

KEGG_result <- KEGG_diff2@result
#以富集结果表Top20为例：
KEGG_top20 <- KEGG_result[1:20,]

#指定绘图顺序（转换为因子）：
KEGG_top20$pathway <- factor(KEGG_top20$Description, levels = rev(KEGG_top20$Description))

#Top20富集数目条形图：
mytheme <- theme(axis.title = element_text(size = 13),
                 axis.text = element_text(size = 11), 
                 plot.title = element_text(size = 14, hjust = 0.5, face = "bold"), 
                 legend.title = element_text(size = 13), 
                 legend.text = element_text(size = 11)) #自定义主题

p <- ggplot(data = KEGG_top20, 
            aes(x = Count, y = pathway, fill = -log10(pvalue)))+
  geom_bar(stat = "identity", width = 0.8) + 
  scale_fill_distiller(palette = "RdYlBu",direction = 1) +
  labs(x = "Number of Gene",
       y = "pathway",
       title = "KEGG enrichment barplot") +
  theme_bw() +
  mytheme
p

p <- barplot(go, color = "p.adjust", showCategory = 20,split="ONTOLOGY") +facet_grid(ONTOLOGY~., scale="free")


GO_result <- go@result
#以富集结果表Top20为例：
GO_top20 <- GO_result[1:10,]

#指定绘图顺序（转换为因子）：
GO_top20$pathway <- factor(GO_top20$Description, levels = rev(GO_top20$Description))


ggplot(data = GO_top20, 
            aes(x = Count, y = pathway, fill = -log10(pvalue)))+
  geom_bar(stat = "identity", width = 0.8) + 
  scale_fill_distiller(palette = "RdYlBu",direction = 1) +
  labs(x = "Number of Gene",
       y = "pathway",
       title = "GO enrichment barplot") +
  theme_bw() +
  mytheme



library(ReactomePA)
## ReactomePA v1.42.0  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
## 
## If you use ReactomePA in published research, please cite:
## Guangchuang Yu, Qing-Yu He. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Molecular BioSystems 2016, 12(2):477-479

x <- enrichPathway(gene=diff_entrez$ENTREZID, 
                   pvalueCutoff = 0.05, 
                   readable=TRUE)




Rectome_result <- x@result
#以富集结果表Top20为例：
Rectome_top20 <- Rectome_result[1:10,]

#指定绘图顺序（转换为因子）：
Rectome_top20$pathway <- factor(Rectome_top20$Description, levels = rev(Rectome_top20$Description))


ggplot(data = Rectome_top20, 
       aes(x = Count, y = pathway, fill = -log10(pvalue)))+
  geom_bar(stat = "identity", width = 0.8) + 
  scale_fill_distiller(palette = "RdYlBu",direction = 1) +
  labs(x = "Number of Gene",
       y = "pathway",
       title = "Reactome enrichment barplot") +
  theme_bw() +
  mytheme
